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We present some new results regarding simulations of finite density QCD based on a canonical 
approach. A previous study has shown that such simulations are feasible, at least on small lattices. 
In the current study, we investigate some of the issues left open: we study the errors introduced by 
our approximation of the action and we show how to tune it to reduce the cost of the simulations 
while keeping the errors under control. To further reduce the cost of the simulations, we check 
the reliability of reweighting method with respect to the baryon number. Finally, using these 
optimizations, we carry out the simulations at larger densities than in our previous study to look 
for signals of a phase transition. 
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1. Introduction 

In recent years, full QCD simulations have become feasible due to development of new al- 
gorithms and increasing computational power. Lattice simulations using dynamical fermions can 
now be performed at finite temperature and zero baryon density. However, simulations at non-zero 
density remain a challenge for lattice QCD due to the complex nature of the fermionic determinant 
where the conventional Monte Carlo methods fail. The standard solution of splitting the action into 
a real an positive part and a phase fails due to the so called sign problem and overlap problem. 
A method to address the overlap problem has been proposed This method uses the canonical 
partition function as a starting point, in contrast to the more common approach based on the grand 
canonical partition function. Using this method, simulations that employ an exact calculation of 
the determinant [Q] were carried out. It was shown that finite density simulations based on the 
canonical partition function are feasible and a program was outlined to look for the phase transition 
at non-zero baryon densities. New simulations were carried out to address some of the issues not 
covered in the previous study; in this paper we report the results of these simulations. 

Due to practical considerations, simulations based on the canonical partition function have to 
use an approximation of the action. The original study looked at deviations from the expected quark 
number in the box to gauge the errors due to this approximation. In this work, we use a more direct 
method: we measure the same physical quantities while reducing the errors of our approximation. 
We find that the errors are small except for the chemical potential measured at the symmetric point 
of discrete Fourier transform. As it turns out, the source of these errors is very easy to understand 
and correct. This also leads us to a very simple criterion that can be used to reduce the simulation 
cost while keeping the errors at a minimum. 

To further reduce the cost of our simulations we employ reweighting. We analyze the reliability 
of reweighting in baryon number using the methods we tested in [Qj. By comparing the results of 
reweighting with the direct simulations we find that we can extrapolate reliably by at least one 
baryon. 

Finally, we use the ensembles generated for these studies to look for a signal for phase transi- 
tion at large baryon number. We measure both the chemical potential and the plaquette distribution 
and vary the baryon number while keeping the temperature fixed. We scan the density space at two 
different temperatures close to the critical temperature. We don't find any clear signal for a phase 
transition on the lattices we studied. 

2. Canonical approach 

Finite density ensembles were generated using the canonical approach To build canonical 
partition function, we start from the fugacity expansion of the grand canonical partition function 

Z(V,T,ll)=Y,Zc(V,T,k)e^ T , (2.1) 

k 

where k is the net quark number and Zq is the canonical partition function of the system. On the 
lattice, we can easily compute the Fourier transform of the grand canonical partition function 

Z(V,T,n) = J &\J&\ff&\if^ Ss(uySf(p '- U ^' v) (2.2) 
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with imaginary chemical potential to get the canonical partition function 

1 f 2n 

Z c (V,T,k) = —J o d((,e-^Z(V,T^)\^ T . (2.3) 

We will specialize this to the case of two degenerate flavors. After integrating out the fermionic 
part, we get a simple expression 

Z c (V,T,k)= f ®Ue- s s iu) det k M 2 (U), (2.4) 

where 

1 r 2n 

det k M 2 (U) = —J o d§e- lk *feM(m,ii;U)\= m , (2.5) 

is the projected determinant with the fixed net quark number k. Zq is the starting point for our 
simulations. For practical reasons, we replace the continuous Fourier transition with its discrete 
version 

tel k M\U) = 1 N fe-^detM(U^ 2 , 0; = (2.6) 

N 7=0 N 

The partition function that we use in our simulations is: 

Z c (V,T,k) = J DUe^ {u) Redet k M 2 (U). (2.7) 

For more details, we refer the reader to the original paper |0|. 

We want to emphasize that the determinant of the fermionic matrix needs to be calculated 
in every step while generating canonical ensembles It is well known that the calculation of 
determinant calculation is very expensive even for a small lattice (on a 4 4 lattice, the dimension 
of the matrix is 3072). It is also possible to simulate this action using a noisy estimator for the 
determinant (Noisy Monte Carlo [||]). In this paper, we show results based only on the exact 
determinant calculation. 



3. Discretization errors and optimized Fourier transform 

To compute the Fourier transform of the determinant, we have to evaluate the determinant for 
all possible phases. This is clearly not feasible and we are forced to use a discrete version of the 
Fourier transform. As shown in the errors introduced by this approximation are proportional 
to exp(— [F(\N\ —n) —F{n)]/T). It is easy to see that as N increases the contamination decreases 
exponentially, while the cost increases only linearly. We can then reduce the error for a modest 
increase in computation time. The only problem is to estimate the effect of the approximation. 
In [Q], an indirect method was proposed that gauged the effect of contamination by measuring the 
difference between the expected number of baryons in the box and the measured one. 

In this study we use a direct approach: we increase ,/V and measure the observables of interest 
again. One of them is the absolute value of the Polyakov loop: 

,| f u (i f l")o (31) 

m) (o> ' {5A) 
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where 



a(U) 



Redet k M 2 (U) 



Redet k M 2 {U) 



(3.2) 



is the phase and () stands for the average over the ensemble generated with measure \Redet k M 2 (U) 
The other one is the baryon chemical potential 



F(n B + l)-F(n B ) 



("B + l 



■n B 



where 



r(u) 



1 Zc(3n B + 3) 
J3 n Z c {3n B ) 



Redet 3 „ s+3 M 2 ([/) 



i . (ygO), 



Redet 3;ifi M 2 (£/) 



(3.3) 



(3.4) 



In the equations above, N appears implicitly in the definition of the probability measure for 
the generated ensemble and in the definition of the observables. The reason we use the same value 
in these two definitions is that we want to have the generated ensemble as close as possible to the 
target one. However, it is possible to do reweighting by measuring the observables on a target 
ensemble with N = N a , on ensembles generated with N = N e . In other words, N e is used in defining 



the measure 



Redet^M" 



and in the definition of the denominators for a(U) and j(U). N is used 
in defining the numerators for a(U) and y(U). The expressions for the Polyakov loop and the 
chemical potential remain the same in terms of Oc(U) and y(U). 

The reason for this separation is that it is expensive to increase N e since we have to evaluate 
the determinant N e times at every accept-reject step; whereas, N is relatively cheaper to increase 
since we measure the projected determinant based on N only on the saved configurations (a factor 
of 100 to 1000 times less evaluations than in the previous case). On the other hand, we do not want 
to have N much larger than N e since this can lead to an overlap problem. 

To investigate the effects of our approximation we first increased N from 12, the value used in 
our previous study [0|, to 24. In this case, we could directly use the ensembles already generated. 
In Table [I], we present the results of this comparison for j3 = 5.20 on a 4 4 lattice with m % rj lGeV. 
We see that the Polyakov loop measurements are not affected by this change, while the chemical 
potential measurements at k = 3 disagree. To understand this result we decided to carry out a more 
direct check by generating an ensemble with N e = 24. 





Polyakov loop 


Baryon chemical potential (n nB /T) 


k 


N e>0 = 12 


N e = 12N = 24 


N e , = 24 


N e . = 12 


N e = \2N = 2A 


N e>0 = 24 





0.192(3) 


0.192(3) 


0.193(8) 


3.40(21) 


3.40(21) 


3.20(15) 


3 


0.363(6) 


0.363(6) 


0.354(13) 


3.09(11) 


3.78(11) 


3.74(10) 



Table 1: Polyakov loop and baryon chemical potential: direct measurements vs reweighting. 

In Table [j] the results from the direct simulations are presented. We see that the Polyakov 
loop is the same for all three cases. The baryon chemical potential for k = 3 is very interesting: it 
agrees with the value for N e = 12 and Af c = 24 rather than the one produced with N e = N = 12. To 
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understand this, we first note that 

Z c (V,T,k)= £ Zc(V,T,k + mN), (3.5) 

m=—°° 

where Zq (V, T, k + mN) is defined using the continuous Fourier transform. Using the above property 
together with the fact that Zc(V,T,k) = Zc(V,T,—k), we get for the symmetric point (k = 6 for 
# = 12) 

Z C (V,T,6) = £ Z c (V,T,6 + m\2) 

m=— °° 

= ... +Z C (V,T, -6) +Z C (V,T, 6) +Z C (V,T, 18) + ... 

« Z c (V,r,-6)+Z c (V,r,6) =2Zc(V,r,6). (3.6) 



Using Eq. ( 3.3) we get 



WES 1 / 1, ^ ->°g^f - = /r "'" 2 <37) 

The ln2 correction changes 3.09(11) inTable[I]to3.78(ll) which agrees with the result ofreweight- 
ing and direct simulation. 

Although the baryon chemical potential measurement at the symmetric point could be cor- 
rected by an extra In 2, the contamination from the other sector can introduce errors in other ob- 
servables. To be on the safe side it is better to carry out simulations with N e larger than twice the 
value of k. Since the contamination errors seem much smaller than the statistical error we will carry 
out our simulations with N e = 2k + 3. Before we move on, we want to point out that these conclu- 
sions might change when we lower the mass of the quark as pointed out in [Q] the contamination is 
suppressed by a factor proportional to the baryon free energy which is expected to decrease as we 
lower the quark mass. 



4. Reliable range of reweighting 

In the previous section, we have shown that we can generate ensembles based on a smaller 
N e and use a form of reweighting to get results corresponding to a different ensemble. We can 
generalize this reweighting to change not only N to N but also k from the value used to generate 
the ensemble, k e , to the one that will be used in the observable, k . The only difference is that 
one needs to adjust the phase factors accordingly. To determine the reliability of the method, we 
compare the results of reweighting to the ones derived from the direct measurements. In Fig. [j], we 
present the results based on 3 ensembles with k e = 0, 3 and 6. They are all generated at j3 = 5.20, 
using N e = 12 and the same quark mass. We use the reweighting to produce results corresponding 
to N = 48 and k = 0, 3, 6, 15 and we compare these results with those of direct measurements. 
For the direct measurements we use k e = k and N e = N and we set N e = 24 for k e = 0, 3 and 6, 
N e = 21 for k e = 9 and N e = 39 for k e = 12 and 15. 

We find that reweighting is reliable for quite a large range of baryon numbers. The fact that we 
can extrapolate from k e = to k = 9 is quite surprising but we cannot thrust this property to remain 
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Polyakov loop Baryon chemical potential 




Figure 1: Left: Polyakov loop from reweighting against direct simulations. Right: The same for baryon 
chemical potential. k () = k e except for the case of reweighting. 

true at different temperatures and volumes. However, our plan was to use reweighting to interpolate 
rather than extrapolate: we plan to scan the baryon number space sparsely (skipping some values 
of k) and filling these gaps by reweighting. This study suggests that this method should produce 
reliable results if \k e — k \ is no larger than 6. 

5. Possible phase diagram 

We now turn to the physical problem - the QCD phase diagram. It is expected that two flavor 
QCD phase diagram will have a critical point at non-zero chemical potential where the first order 
phase transition turns to a crossover. In the canonical ensemble, the first order phase transition 
will be represented by a coexistence region. To search for the phase boundaries we scan the phase 
diagram by varying the density while keeping the temperature fixed. The baryon chemical potential 
should exhibit an "S-shape" [^J as one crosses the coexistence region. The Maxwell construction 
can then be used to determine the phase boundaries. Another possible signal can be detected by 
looking at the histogram of the plaquette. ^ 

To test these ideas, we have generated two sets of ensembles at two temperatures close to the 
critical temperature. In Fig. ||, we plot the baryon chemical potential and the histogram of the 
plaquette. We notice a peak shift in the histogram plot but the chemical potential does not exhibit 
a clear "S-shape". Although there is no evidence of an "S-shape", the peak shift suggests that 
a transition happens between 3.6 to 18 times the normal nuclear matter density whereas Ph. de 
Forcrand and S. Kratochvila found the range to be 1 to 10 times of the nuclear density The 
difference may be due to the heavy quark mass we use in our study and the fact that they studied a 4 
flavor QCD. Our heavy quark mass may also have suppressed the fluctuations and hence obscured 
the "S-shape". 

6. Summary 

In this paper, we present new results that address some of the issues not covered in our previous 
study. We show that the errors introduced by using a discrete Fourier transform are smaller that 
the statistical errors with one exception that we analyze and show how to correct. We find that 
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Figure 2: Chemical potential and plaquette histograms for temperatures close to the transition temperature. 

the optimal strategy is to set N e = 2k + 3 to reduce the simulation cost while keeping the errors 
under control. We also analyze the reliability of reweighting in baryon number by comparing the 
results of reweighting with direct simulations. We find that the method is reliable for a larger range 
than we expected and this will allow us to further reduce the cost of our simulations in the future. 
Finally, we have used the ensembles generated for this study to look for the phase transition at 
larger baryon numbers. Although the chemical potential plot does not show a clear signal for phase 
transition, the peak shift in the plaquette histogram from small densities to large densities suggests 
a transition. 

We plan to continue this study with simulations on a larger lattice (6 3 x 4) where we hope to 
be able to detect a phase transition. We also plan to use improved gauge and fermionic actions to 
run simulations at smaller quark masses. 
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